library(tidyverse)
library(colorDF)
library(ggplot2)
library(ggbeeswarm)
library(tmod)
library(colorDF)
library(cowplot)
theme_set(theme_minimal())
This is a data set that I used in my poster. It is available from GEO and contains RNA-Seq counts for healthy individuals and COVID-19 patients, described by Mick et al. (2020).
The data set has been already processed and downloaded; however, below you will find the code to replicate this procedure (if you want).
## this sometimes fails. The RDS of the object is provided
Sys.setenv(VROOM_CONNECTION_SIZE=8*131072)
## with getGEO, we can only get the phenoData
geo <- getGEO("GSE156063")[[1]]
saveRDS(geo, file="GSE156063.geo")
geo <- readRDS("GSE156063.geo")
covar <- as(phenoData(geo), "data.frame")
## remove columns with only one value
boring <- map_lgl(covar, ~ length(unique(.x)) == 1)
covar <- covar[ , !boring ]
## clean up
covar <- covar %>%
dplyr::rename(gender = "gender:ch1") %>%
mutate(disease = gsub(" .*", "", .data[["disease state:ch1"]])) %>%
mutate(label = description) %>%
mutate(group = disease) %>%
arrange(description) %>%
dplyr::select(all_of(c("title", "label", "gender", "disease", "group")))
## the counts must be downloaded from GEO separately.
if(!file.exists("GSE156063_swab_gene_counts.csv.gz")) {
download.file("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE156nnn/GSE156063/suppl/GSE156063_swab_gene_counts.csv.gz",
"GSE156063_swab_gene_counts.csv.gz")
}
counts <- read_csv("GSE156063_swab_gene_counts.csv.gz")
.tmp <- counts[[1]]
counts <- as.matrix(counts[,-1])
rownames(counts) <- .tmp
#counts <- counts[ , covar$description ]
counts <- counts[ , covar$label ]
lcpm <- edgeR::cpm(counts, log=TRUE)
#stopifnot(all(colnames(counts) == covar$description))
stopifnot(all(colnames(counts) == covar$label))
annot <- data.frame(ENSEMBL = rownames(counts))
.tmp <- mapIds(org.Hs.eg.db, annot$ENSEMBL, column=c("ENTREZID"), keytype="ENSEMBL")
annot$ENTREZID <- .tmp[ match(annot$ENSEMBL, names(.tmp)) ]
.tmp <- mapIds(org.Hs.eg.db, annot$ENSEMBL, column=c("SYMBOL"), keytype="ENSEMBL")
annot$SYMBOL <- .tmp[ match(annot$ENSEMBL, names(.tmp)) ]
.tmp <- mapIds(org.Hs.eg.db, annot$ENSEMBL, column=c("GENENAME"), keytype="ENSEMBL")
annot$GENENAME <- .tmp[ match(annot$ENSEMBL, names(.tmp)) ]
sel <- covar$group %in% c("no", "SC2")
counts <- counts[ , sel ]
covar <- covar[ sel, ]
covar$group <- as.character(covar$group)
covar$group.disease <- paste0(covar$group, '_', covar$disease)
rownames(covar) <- covar$label
saveRDS(covar, file="covar.rds")
saveRDS(counts, file="counts.rds")
## DESeq2 calculations
## manual cache, since the operation takes a long time
sel <- !is.na(covar$group)
ds2 <- DESeqDataSetFromMatrix(counts[,sel], colData=covar[sel, ], design=~ disease)
ds2_file <- "ds2_cached.rds"
if(!file.exists(ds2_file)) {
message("Running DESeq2")
ds2 <- DESeq(ds2)
saveRDS(ds2, file=ds2_file)
} else {
message("Reading ds2 from manual cache")
ds2 <- readRDS(ds2_file)
}
First, take a look at the structure of the data:
covar <- readRDS("data/covar.rds") # covariates
summary_colorDF(covar)
The results table is stored under “data/disease_SC2_vs_no.tsv”. In
data/annotation.tsv, there is the annotation of the genes
in this data set.
dsres <- read_table("data/disease_SC2_vs_no.tsv")
annot <- read_table("data/annotation.tsv")
dsres <- merge(annot, dsres, by="ENSEMBL")
Basic CERNO analysis:
gl <- dsres$SYMBOL[ order(dsres$pvalue) ]
tres <- tmodCERNOtest(gl)
print(tres)
## # tmod report (class tmodReport) 8 x 107:
## # (Showing rows 1 - 5 out of 107)
## │ID │Title │cerno│N1 │AUC │cES │P.Value
## DC.M3.4│ DC.M3.4│Interferon │ 425│ 48│ 0.89│ 4.4│<2e-16
## DC.M1.2│ DC.M1.2│Interferon │ 296│ 24│ 0.99│ 6.2│<2e-16
## DC.M5.12│DC.M5.12│Interferon │ 356│ 56│ 0.82│ 3.2│<2e-16
## DC.M3.2│ DC.M3.2│Inflammation │ 528│ 112│ 0.82│ 2.4│<2e-16
## LI.M75│ LI.M75│antiviral IFN signature│ 198│ 21│ 0.91│ 4.7│<2e-16
## │adj.P.Val
## DC.M3.4│ 3.9e-40
## DC.M1.2│ 5.6e-35
## DC.M5.12│ 7.7e-25
## DC.M3.2│ 1.8e-24
## LI.M75│ 5.6e-20
Which we can quickly visualize:
ggPanelplot(tres)
We can also add information about whether the genes go up or down.
Note that we need to put tres inside a list, so
ggPanelplot can match tmod results with the
sgenes object.
sgenes <- tmodDecideTests(dsres$SYMBOL, lfc = dsres$log2FoldChange,
pval=dsres$pvalue)
names(sgenes) <- "SC2_vs_no"
tres <- list("SC2_vs_no"=tres)
ggPanelplot(tres, sgenes=sgenes)
Another built-in gene set are the cell surface markers.
data(cell_signatures)
tres2 <- tmodCERNOtest(gl, mset=cell_signatures)
tres2 <- list("SC2_vs_no"=tres2)
sgenes2 <- tmodDecideTests(dsres$SYMBOL, lfc = dsres$log2FoldChange,
pval=dsres$pvalue, mset=cell_signatures)
names(sgenes2) <- "SC2_vs_no"
ggPanelplot(tres2, sgenes=sgenes2, mset=cell_signatures)
The mset parameter is used in most of the functions to use another gene sets.
Let us consider the gene set “LI.M68” – RIG-1 like receptor signalling. First, show the genes in this gene set:
data(tmod) # make the built-in gene sets visible
getModuleMembers("LI.M68", mset=tmod)
## $LI.M68
## [1] "DDX58" "IRF7" "DHX58" "TRIM25" "NFKBIA" "IFIH1" "TNF" "ISG15"
## [9] "CXCL10" "IL8"
Evidence plot – a ROC curve:
par(mfrow=c(2,2))
evidencePlot(gl, m = "DC.M1.2", gene.labels=TRUE)
evidencePlot(gl, m = "LI.M68", gene.labels=TRUE)
evidencePlot(gl, m = "LI.M165")
evidencePlot(gl, m = "LI.M4.0")
The basic idea is as follos: run PCA only on a selected group of
genes, use PC1 instead of average gene expression of that gene set (mind
the signs!). Below, counts.rds contains the raw RNA-seq
counts, which I first convert to log2 counts per million.
counts <- readRDS("data/counts.rds")
lcpm <- edgeR::cpm(counts, log=TRUE)
mm <- getModuleMembers("LI.M68", mset=tmod)[[1]]
sel <- annot$SYMBOL %in% mm
covar$eig <- prcomp(t(lcpm[ sel, ]), scale.=TRUE)$x[,1]
We can use eigengenes to nicely plot the results:
ggplot(covar, aes(x=disease, y=eig)) + geom_boxplot(outlier.shape = NA) +
geom_beeswarm(cex=3)
(However, mind the sign!)
In tmod, eigengenes can be computed directly. Tmod makes sure that the eigengene correlates positively with the majority of the genes in the gene set and reverses the sign when necessary:
selm <- tres[[1]]$ID[1:4]
eig <- eigengene(lcpm, annot$SYMBOL)
eig <- as.data.frame(t(eig)) %>%
select(all_of(selm)) %>%
rownames_to_column("label") %>%
inner_join(covar, by="label")
plots <- map(selm, ~ ggplot(eig, aes_string(x="disease", y=.x)) +
geom_boxplot(outlier.shape=NA) +
geom_beeswarm(cex=3) +
ggtitle(tres[[1]]$Title[ match(.x, tres[[1]]$ID) ]))
plot_grid(plotlist=plots)
Using a method like tmodCERNOtest, tmodUtest or fsgsea makes it possible to apply gene set enrichment to any single ranked list of genes.
(I have selected the components based on correlation with
covar$disease)
pca <- prcomp(t(lcpm), scale.=TRUE)
pca_df <- pca$x %>% as.data.frame %>%
rownames_to_column("label") %>% merge(covar)
ggplot(pca_df, aes(x=PC5, y=PC6, color=disease)) + geom_point()
We can sort all the genes by their decreasing absolute loadings in the pca object.
pcares <- map(1:10, ~ {
gl <- annot$SYMBOL[ order(-abs(pca$rotation[,.x])) ]
tmodCERNOtest(gl)
})
names(pcares) <- colnames(pca$rotation)[1:10]
sgenes <- tmodDecideTests(annot$SYMBOL,
lfc=pca$rotation[, 1:10],
pval=NULL,
lfc.thr = 0)
ggPanelplot(pcares, sgenes=sgenes, filter_row_auc = .85)
ggplot(pca_df, aes(x=PC2, y=PC3, color=disease)) + geom_point()
Check the function tmod_pca for more details.
Converting the data frame from the msigdbr package for use with tmod is straightforward:
library(msigdbr)
msig_df <- msigdbr("Homo sapiens")
goset <- msig_df %>%
filter(gs_subcat == "GO:BP") %>%
makeTmodFromDataFrame(module_col="gs_name",
title_col="gs_name",
feature_col="gene_symbol")
library(DT)
tmodCERNOtest(gl, mset=goset, qval = 1e-3) %>%
filter(N1 > 5) %>%
datatable() %>%
formatSignif(c("adj.P.Val", "P.Value")) %>%
formatRound(c("cES", "AUC"))
tres[[1]] %>% filter(N1 < 20) %>% pull(ID) %>% upset()
foo <- tres[[1]] %>%
slice(1:20) %>%
pull(ID) %>%
modCorPlot(stat = "overlap")
foo
Which genes are driving the enrichment?
evidencePlot(gl, m = "LI.M165", mset=tmod)
evidencePlot(gl, m = "LI.M165", mset=tmod, style="gsea")
Get the leading edge:
lea <- tmodLEA(gl, "LI.M165")
lea
## $LI.M165
## [1] "IFIT1" "IFI27" "CXCL11" "RSAD2" "SIGLEC1" "PARP9"
## [7] "LAMP3" "IFIT2" "IFIT3" "IFIH1" "HERC5" "CCL8"
## [13] "SERPING1" "HESX1" "CCL20" "DSE" "IFNGR2" "RBM47"
## [19] "MX2" "BCL2A1" "TNFAIP6" "FAM20C" "CSF2RB" "DOCK4"
## attr(,"N")
## [1] 35
## attr(,"LEA")
## [1] 2831
## attr(,"LEA.frac")
## [1] 0.6857143
##
## attr(,"class")
## [1] "tmodLEA" "list"
clr <- ifelse(gl %in% lea[[1]], "red", "grey")
names(clr) <- gl
evidencePlot(gl, m = "LI.M165", mset=tmod,
gene.labels=TRUE,
gene.colors = clr)
And GSEA style:
clr <- ifelse(gl %in% lea[[1]], "red", "grey")
names(clr) <- gl
evidencePlot(gl, m = "LI.M165", mset=tmod,
gene.labels=TRUE,
gene.colors = clr,
style="gsea")
abline(v=attr(lea[[1]], "LEA"), col="red", lwd=2)
This example data set has been published by Weiner et al. (2019). RNA has been collected daily from individuals vaccinated with one of seven vaccines over the course of several days. The data set has been collected using microarrays. The pre-calculated results include two influenza vaccines, one containing an adjuvant and the other not. For details how to get this table, see the online tmod manual here, section “Transcriptional responses to vaccination”.
tt <- readRDS("data/vaccination_toptable.rds")
The subset contains results for 5 time points for two vaccines (10
contrasts in total). The two vaccines are denoted with “A” (no adjuvant)
and “F” (adjuvant). We will run gene set enrichment on each of these
time points using tmodCERNOtest.
contrasts_v <- paste0(rep(c("F", "A"), each=5), "_", 1:5)
sort_cols <- paste0("qval.", contrasts_v)
res <- map(sort_cols, ~ {
tmodCERNOtest(tt$GENE_SYMBOL[ order(tt[[.x]]) ])
})
names(res) <- contrasts_v
The figure below shows the results.
ggPanelplot(res, filter_row_q = 1e-3, filter_row_auc = .85, q_thr=0.01)
Do these gene sets contain genes that go up or down? For this, we need to determine how many genes in each gene set are significantly going up or down.
pval <- tt %>% select(starts_with("qval"))
lfc <- tt %>% select(starts_with("logFC"))
sgenes <- tmodDecideTests(tt$GENE_SYMBOL, lfc=lfc, pval=pval)
names(sgenes) <- contrasts_v
ggPanelplot(res, sgenes=sgenes, filter_row_q = 1e-3, filter_row_auc = .85, q_thr=0.01)
Are responses on day 1 in Fluad similar to responses on day 1 in Agrippal?
\[d_s = \log_2 FC_A \cdot \log_2 FC_B \cdot |\log_{10}p_A + \log_{10}p_B|\]